subroutine set_ic( dfu, uxMat, f0vec, fvc0T, it)    ! outputs

  use dimensions

  use fft_wrap, only   : fft_forward, fft_backward

  use parameters, only : kz, vPhase, u, &                             ! real,     phys
                         nz, iWaveFrame, iModx, iModz, &              ! integer,  phys
                         ModAmp,  &                                   ! real,     phys
                         restart, &                                   ! integer,  num
                         idfMod, &                                    ! integer,  num 
                         pi, iunit, isign, &                          ! constants
                         runname

  implicit none

!-- interface input/output --

  real,    intent(out) :: dfu(jz, jv, jx, jDirection)
  real,    intent(out) :: uxMat(jx, jz, jDirection)
  real,    intent(out) :: f0vec(jv)
  complex, intent(out) :: fvc0T(jv)
  integer, intent(out) :: it


!-- local variables ---

  integer              :: ix, iz, iv, iDirection  ! counters in r-space  
  real                 :: MoDinZ(jz), MoDinX(jx) 
  real                 :: v, f0


  f0(v)=( 1/sqrt(2*pi) )*exp(-.5*(v+vPhase*iWaveFrame)*(v+vPhase*iWaveFrame))

!-- create vector "f0vec" --


  do iv=1,jv
     v = vvec(iv)
     f0vec(iv) = f0(v)   
  end do
                 
  fvc0T(:)=f0vec(:)      
  call fft_forward(jv,fvc0T,fvc0T)   !  fvc=FT(dfu)
  fvc0T=fvc0T/float(jv)             !  normalized velocity  FT of equilibriium


!-- read restart or create new --

  if (restart.ne.0) then         ! read restart info

       write(*,*)
       write(*,*) 'reading restart'

       call read_bin(dfu, uxMat, it)        

  else                   ! start fresh

       write(*,*)
       write(*,*) 'fresh start'

       it = 0
       dfu=0.0

       do iDirection = 1, 2
          uxMat(:,:,iDirection) = u * iSign(iDirection)
       enddo
      
       if(idfMod==1)then

          if( iModx.ne.1 ) iModx=0 ;   if( iModz.ne.1) iModz=0
         
          MoDinZ(:) =    iModz * cos(  kz * zvec(:) )
          MoDinX(:) =    iModx * cos( dkx * xvec(:) )
         
          do iDirection = 1,2
             do iv=1,jv                                                      
                do ix=1,jx
                   do iz=1,jz
                      dfu(iz,iv,ix,iDirection) =  ModAmp * f0vec(iv) *  (  MoDinZ(iz) + MoDinX(ix) )
                   enddo
                enddo
             enddo
          enddo !   iv and iDirection

       endif ! if idfMod


  end if  !  end restart test


end subroutine set_ic
